#### Check packages

requiredPackages <- c("CVXR")
for (p in requiredPackages) {
  if (!require(p, character.only = TRUE)) {
    install.packages(p, dependencies = TRUE)
    library(p, character.only = TRUE)
  }
}



################################################################################
################################################################################

#### Mean-CVaR Optimization 

mean_CVaR_optimization <- function(n, nc, N, alpha, S0, C, rf, t, ST, k, D, r_target) {
  
  w <- Variable(n)
  p <- Variable(n, nc)  
  q <- Variable(1)
  z <- Variable(N)
  Er <- Variable(1)  
  x <- Variable(n)
  r <- Variable(N)  
  
  objective <-  q + 1/(N * (1 - alpha)) * sum(z)
  
  constraints <- list()
  
  for (i in 1:N) {
    constraints <- c(constraints, z[i] >= 0)
    constraints <- c(constraints, z[i] + r[i] + q >= 0 )
  }
  
  for (j in 1:n) {
    constraints <- c(constraints, w[j] >= 0)
    for (l in 1:nc) {
      constraints <- c(constraints, p[j,l] >= 0)
    }
  }
  
  
  constraints <- c(constraints, sum(w) == 1)
  
  for (j in 1:n) {
    constraints <- c(constraints, x[j] == w[j] / S0[j] )  
    constraints <- c(constraints, sum(p[j, 1:nc]) - x[j] <= 0)
  }
  

  for (i in 1:N) {
    asset_returns <- 0
    for (j in 1:n) {
      asset_returns <- asset_returns + sum(x[j] * (ST[i,j] + D[j]))
    }
    
    option_returns <- 0
    for (j in 1:n) {
      for (l in 1:nc) {
        option_payoff <- pmax(ST[i,j] - k[j, l], 0)
        option_returns <- option_returns + p[j, l] * (C[j, l] * exp(rf * t) - option_payoff)
      }
    }
    
    constraints <- c(constraints, r[i] == asset_returns + option_returns - 1)
  }
  
  constraints <- c(constraints, Er == sum(r) / N)
  constraints <- c(constraints, r_target <= Er)
  
  problem <- Problem(Minimize(objective), constraints)
  result <- solve(problem, solver = "GLPK")
  
  details <- list(
    w = result$getValue(w),
    p = result$getValue(p),
    x = result$getValue(x),
    q = result$getValue(q),
    status = result$status,
    r = result$getValue(r),
    value = result$value
    
  )
  
  return(details)
}


################################################################################
################################################################################

#### Variance Optimization

Variance_optimization <- function(n, nc, N, S0, C, rf, t, ST, k, D, r_target) {
  
  w <- Variable(n)
  p <- Variable(n, nc)  
  q <- Variable(1)
  Er <- Variable(1)  
  x <- Variable(n)
  r <- Variable(N)  
  
  objective <-  1/(N-1)*sum((r - sum(r)/N)^2)
  
  constraints <- list()
  
  for (j in 1:n) {
    constraints <- c(constraints, w[j] >= 0)
    for (l in 1:nc) {
      constraints <- c(constraints, p[j,l] >= 0)
    }
  }
  
  
  constraints <- c(constraints, sum(w) == 1)
  
  for (j in 1:n) {
    constraints <- c(constraints, x[j] == w[j] / S0[j] )  
    constraints <- c(constraints, sum(p[j, 1:nc]) - x[j] <= 0)
  }
  

  for (i in 1:N) {
    asset_returns <- 0
    for (j in 1:n) {
      asset_returns <- asset_returns + sum(x[j] * (ST[i,j] + D[j]))
    }
    
    option_returns <- 0
    for (j in 1:n) {
      for (l in 1:nc) {
        option_payoff <- pmax(ST[i,j] - k[j, l], 0)
        option_returns <- option_returns + p[j, l] * (C[j, l] * exp(rf * t) - option_payoff)
      }
    }
    
    constraints <- c(constraints, r[i] == asset_returns + option_returns - 1)
  }
  
  constraints <- c(constraints, Er == sum(r) / N)
  constraints <- c(constraints, r_target <= Er)
  
  problem <- Problem(Minimize(objective), constraints)
  result <- solve(problem, solver = "ECOS")
  
  details <- list(
    w = result$getValue(w),
    p = result$getValue(p),
    x = result$getValue(x),
    status = result$status,
    r = result$getValue(r),
    value = result$value
    
  )
  
  return(details)
}


################################################################################
################################################################################

#### Semivariance Optimization

Semivariance_optimization <- function(n, nc, N, S0, C, rf, t, ST, k, D, r_target) {
  
  w <- Variable(n)
  p <- Variable(n, nc)  
  q <- Variable(1)
  z <- Variable(N)
  Er <- Variable(1)  
  x <- Variable(n)
  r <- Variable(N)  
  
  objective <-  1/N * sum(z^2)
  
  constraints <- list()
  
  for (i in 1:N) {
    constraints <- c(constraints, z[i] >= 0)
    constraints <- c(constraints, z[i] + r[i] + 1 - exp(rf*t) >= 0 )
  }
  
  for (j in 1:n) {
    constraints <- c(constraints, w[j] >= 0)
    for (l in 1:nc) {
      constraints <- c(constraints, p[j,l] >= 0)
    }
  }
  
  
  constraints <- c(constraints, sum(w) == 1)
  
  for (j in 1:n) {
    constraints <- c(constraints, x[j] == w[j] / S0[j] )  
    constraints <- c(constraints, sum(p[j, 1:nc]) - x[j] <= 0)
  }
  
  for (i in 1:N) {
    asset_returns <- 0
    for (j in 1:n) {
      asset_returns <- asset_returns + sum(x[j] * (ST[i,j] + D[j]))
    }
    
    option_returns <- 0
    for (j in 1:n) {
      for (l in 1:nc) {
        option_payoff <- pmax(ST[i,j] - k[j, l], 0)
        option_returns <- option_returns + p[j, l] * (C[j, l] * exp(rf * t) - option_payoff)
      }
    }
    
    constraints <- c(constraints, r[i] == asset_returns + option_returns - 1)
  }
  
  constraints <- c(constraints, Er == sum(r) / N)
  constraints <- c(constraints, r_target <= Er)
  
  problem <- Problem(Minimize(objective), constraints)
  result <- solve(problem, solver = "ECOS")
  
  details <- list(
    w = result$getValue(w),
    p = result$getValue(p),
    x = result$getValue(x),
    status = result$status,
    r = result$getValue(r),
    value = result$value
    
  )
  
  return(details)
}


################################################################################
################################################################################

#### Generate random data for testing

random_portfolio <- function(n, nc, N, t, method = "GBM") {
  mu <- runif(1, 0.05, 0.15) 
  sigma <- runif(1, 0.1, 0.3) 
  rf <- runif(1, 0.01, 0.05) 
  r_target <- runif(1, 0.04, 0.08) 
  
  S0 <- runif(n, 50, 150)
  
  C <- matrix(runif(n*nc, 1, 5), n, nc)
  
  k <- matrix(runif(n*nc, 80, 120), n, nc)
  
  D <- runif(n, 1, 3)
  
  ST <- matrix(nrow = N, ncol = n)
  
  if (method == "GBM") {
    for (i in 1:N) {
      Z <- matrix(rnorm(n), nrow = 1)
      ST[i, ] <- S0 * exp((mu - 0.5 * sigma^2) * t + sigma * sqrt(t) * Z)
    }
  } else if (method == "Normal") {
    for (i in 1:N) {
      ST[i, ] <- S0 + rnorm(n, mean = mu * t, sd = sigma * sqrt(t))
    }
  } else {
    stop("Invalid method. Choose 'GBM' or 'Normal'.")
  }
  
  return(list(mu = mu, sigma = sigma, rf = rf, r_target = r_target, S0 = S0, C = C, k = k, D = D, ST = ST,n=n,nc=nc,N=N,t=t))
}


################################################################################
################################################################################


#### Improved Mean-CVaR Optimization (Not working !!!!!!)

improved_mean_CVaR_optimization <- function(n, nc, N, alpha, S0, C, rf, t, ST, k, D, r_target) {
  w <- Variable(n)
  p <- Variable(n, nc)  
  q <- Variable(1)
  z <- Variable(N)
  Er <- Variable(1)  
  x <- Variable(n)
  r <- Variable(N)
  
  objective <-  q + 1/(N * (1 - alpha)) * sum(z)
  
  constraints <- list()
  
  constraints[[length(constraints) + 1]] <- z >= 0
  
  constraints[[length(constraints) + 1]] <- z + r + q >= 0
  
  constraints[[length(constraints) + 1]] <- w >= 0
  
  constraints[[length(constraints) + 1]] <- p >= 0
  
  constraints[[length(constraints) + 1]] <- sum(w) == 1
  
  constraints[[length(constraints) + 1]] <- x == w / S0
  
  for (j in 1:n) {
    constraints[[length(constraints) + 1]] <- sum(p[j,]) - x[j] <= 0
  }
  
  for (i in 1:N) {
    asset_returns <- sum(x * (ST[i,] + D))
    option_returns <- 0
    for (j in 1:n) {
      for (l in 1:nc) {
        option_payoff <- pmax(ST[i,j] - k[j,l], 0)
        option_returns <- option_returns + p[j,l] * (C[j,l] * exp(rf * t) - option_payoff)
      }
    }
    constraints[[length(constraints) + 1]] <- r[i] == asset_returns + option_returns - 1
  }
  
  constraints[[length(constraints) + 1]] <- Er == mean(r)
  
  constraints[[length(constraints) + 1]] <- r_target <= Er
  
  problem <- Problem(Minimize(objective), constraints)
  result <- solve(problem)
  
  details <- list(
    w = result$getValue(w),
    p = result$getValue(p),
    x = result$getValue(x),
    q = result$getValue(q),
    status = result$status,
    r = result$getValue(r),
    value = result$value
  )
  
  return(details)
}


################################################################################
################################################################################

#### Find the largest expected return target

find_r_target <- function(n, nc, ST, C, S0, rf, t, D) {
  CRP <- matrix(nrow = n, ncol = nc)
  V <- vector("numeric", n)  
  M <- vector("numeric", n)  
  
  for (j in 1:n) {
    for (l in 1:nc) {
      CRP[j, l] <- mean(pmax(ST[, j], 0)) - C[j, l] * exp(rf * t)  
    }
  }
  
  for (j in 1:n) {
    V[j] <- min(min(CRP[j, ]), 0)
    M[j] <- (mean(ST[, j]) + D[j] - V[j]) / S0[j]  
  }
  
  r_target <- max(M) - 1
  return(r_target)
}


################################################################################
################################################################################

#### Calculate mean return given w and p

mean_return <- function(x, p, n, N, nc, ST, C, S0, rf, t, D) {
  r <- vector("numeric", N)  
  for (i in 1:N) {
    asset_returns <- sum(x * (ST[i,] + D))
    option_returns <- 0
    for (j in 1:n) {
      for (l in 1:nc) {
        option_payoff <- pmax(ST[i,j] - k[j,l], 0)
        option_returns <- option_returns + p[j,l] * (C[j,l] * exp(rf * t) - option_payoff)
      }
    }
    r[i] = asset_returns + option_returns - 1
  }
  Er = mean(r)
  return(Er)
}

